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Abstract 

We set up a general framework to study the bifurcations of nearly degenerate modes 
of a Maxwell-Bloch laser model, and we apply it specifically to the study of inter- 
actions of two modes of a laser with broken circular symmetry. We use an explicit 
centre manifold reduction to analyse the behaviour of this laser. Complex bifur- 
cation sequences involving single mode solutions, mode locking and mode beating 
regimes are predicted. These sequences are organized by a Hopf bifurcation with 
1 : 1 resonance and a ^-symmetric Takens-Bogdanov bifurcation. Numerical sim- 
ulations of the original system show good agreement with theory. 
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1 Introduction 



Symmetries pose considerable constraints on the configuration of stationary solutions and 
on the time evolution of dynamical systems as well as on their bifurcation behaviour 
[1-4], to the extent that the presence of symmetries can often dictate the existence of 
certain solution branches independently of the specific model under consideration. These 
ideas have been fruitfully used in optics either to classify the patterns produced by ex- 
periments [5] or by models [6-8], and to infer from these the symmetry of the systems 
under investigation. They have also been used to study in great detail the structure of the 
bifurcations of normal form equations with the same symmetry as the experiments [9-11]. 
Despite the considerable success of these approaches, there are several limitations. In the 
first case, the comparison between theory and experiments is qualitative and the bifur- 
cation diagram is not fully explored. In the second case, there is no relation between the 
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parameters used in the normal form analysis and the physical parameters and, further- 
more, it is impossible to find where the normal forms cease to be a good approximation of 
the system studied. A popular alternative is to ignore symmetry altogether and to study 
models, usually nonlinear partial or integro-differential equations, that are derived from 
first principles. An advantage of this approach is that such models can be valid both close 
to and far from the first instability. A disadvantage is that most of the analysis has to 
be performed numerically, which means in practice that the bifurcation diagram is rarely 
analysed completely. 

In this paper we derive normal form equations for a laser model in the presence of a 
typical breaking of cylindrical symmetry of the cavity. We then study in detail the case 
in which the symmetry is broken by astigmatism. In this way we are able to reduce a 
partial differential model to a system of ordinary differential equations and to take full 
advantage of the symmetry in the study of the bifurcation diagram without losing the 
physical meaning of the control parameters. Furthermore, we are able to find the range 
of the control parameters in which the normal form is in quantitative agreement with 
the full model and where the agreement is only qualitative. In the example discussed in 
this paper, we show that the normal forms predict the correct bifurcation structure and 
amplitude of the solutions for small values of the control parameter; for larger values only 
the nature of the bifurcation points is determined (see Figure 9). The bifurcations turn 
out to be organised by a Z 2 -symmetric Takens-Bogdanov point and by a resonant Hopf 
bifurcation. A Takens-Bogdanov bifurcation commonly occurs in this type of problem 
(see [12,13], for example) and the bifurcation diagrams in [14] sketch some key features 
of associated branching behaviour, but here we provide a complete picture with the full 
state diagram (bifurcation set with associated phase portraits) for the normal form of our 
equations for a wide range of physically significant parameters. 

The rest of this paper is organised as follows. In the next section we introduce a specific 
laser model, namely the Maxwell-Bloch two- level laser [15], and carry out an explicit 
centre manifold reduction at the first bifurcation from the trivial solution as the pump 
parameter is increased. We then in Section 3 look more closely at the particular case 
of astigmatic symmetry-breaking from geometric 0(2) symmetry to D 2 symmetry. The 
additional S* 1 phase-shift yields a bifurcation problem with broken 0(2) x S* 1 symmetry, 
much studied in the literature [3,14,16,17]. We find, however, that the thorough analysis 
of responses to linear symmetry-breaking perturbations given in [14] , the treatment most 
relevant to our problem, does not cover all important cases as it assumes the equivalent 
of nonzero detuning (see Section 3) that may not hold in this problem. For the zero 
detuning case we use the normal form for Hopf bifurcation with 1 : 1 resonance [18], and 
note here how a nonlinear coordinate change converts the problem again into one with 
0(2) x S 1 symmetry, but now with the roles of .SO (2), the group of spatial rotation, 
and S 1 interchanged. This has interesting implications for the interpretation of periodic 
solutions in terms of standing or travelling waves. 

Finally, we compare the results with numerical simulations for the original problem. 
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2 Laser normal forms in the presence of generic symmetry breaking 



We assume that the laser under study is a ring cavity gas laser with approximate cylin- 
drical symmetry as shown schematically in Figure 8. This is roughly the same type of 
laser as used in the experiments of Reference [5] and it has also been used in many other 
experimental [19-26] and theoretical [10,27-30] studies of pattern formation in lasers. We 
use here a standard model for such a system [28], incorporating a number of assumptions. 

First we assume that the light is linearly polarised and that the active medium is isotropic, 
so that the electric field amplitude F can be represented by a scalar field. 

Secondly we assume that the active medium can be represented as a collection of two-level 
atoms (Maxwell-Bloch two-level laser), i.e. that only two atomic levels are involved in the 
lasing process: this allows us to use only two functions to represent the atomic medium, 
namely the polarisation and the population inversion. The polarisation P represents the 
induced polarisation field of the atoms, while the population inversion N measures the 
difference in population between the excited and the ground states. 

Thirdly we assume that either the laser is producing a continuous wave of light or, if it is 
pulsed, that the pulses are long compared to the light period. In this case we can write 
the polarisation and the electric field as the product of a rapidly varying phase term and 
a slowly varying amplitude: 



where uja = cIca is the frequency of the atomic transition, with c the speed of light and 
k,A the corresponding wavenumber. The amplitudes F and P are assumed to be varying 
much more slowly than uja (slowly varying approximation) so that the second derivatives 
of F and P with respect to z and t in Maxwell's wave equation can be neglected with 
respect to their first derivatives [31]. 

Finally, we assume that the gain and losses per cavity round trip are small. This allows 
us to neglect the variation of the electric field amplitude F and of the atomic variables 
P and N along the cavity axis, so that the problem is reduced to two spatial dimensions 
{i.e. the centre plane of the active medium) and time. 

Under these hypotheses, we can write the equations for the model as [28] 



F(x, y, z,t)=\ [F(x, y, z, t)e l{kAZ ~^ + c.c." 
P(x, y, z,t)=l \P(x, y, z, t)e i(kAZ - WAt) + c.c.1 , 




(1) 



(2) 



dF 

~dt 
dP 

~dt 



CF + P, 



-P + X F + FN, 



(3) 



(4) 




(5) 
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where the over-bar symbol indicates complex conjugate. The linear operator C is usually 
a partial differential operator which depends upon the laser cavity and it is related to the 
propagation in the empty cavity, i.e. the cavity without the active medium [28]. In these 
equations time t is non-dimensional and is in units of the polarisation decay time. For a 
very large class of cavities, called stable in laser physics, the electric field remains confined 
within the cavity. For these cavities, the operator £ has a system of eigenfunctions A n (x, y) 



CA n (x,y) = (3 n A n (x,y), 



(6) 



which form an orthonormal basis for the space 7i of L 2 complex functions defined in the 
transverse plane. The lasers we consider here are in this class. The functions A n (x,y) are 
called empty cavity modes because they are also eigenfunctions of the operator TZq that 
propagates the field once around the empty cavity from z = back to itself, namely 



TZ c A n (x,y,0) = e^A n (x,y,0). 



(7) 



The complex coefficient f3 n represents the attenuation and phase shift of the mode n 
per unit time. It is related to the eigenvalue of the mode given in (7) by the relation 
Pn = Hn/T c + i5, where T c is the cavity round trip time and 6 is a detuning parameter, 



(8) 



Here uc is a reference frequency with respect to which the propagation phase shifts of 
the modes are measured. The case in which C is an integral operator can be dealt with 
in the same way, but the previous relations are slightly different. 



The energy fed into the laser is represented by the pump parameter x — x{ x -iU) an d 
the decay rate of the population inversion is given by 7. In all that follows we assume 
that x is space- independent. The appropriate calculations for a space-dependent pump 
are sketched in Appendix B. 

Equations (3-5) have a trivial solution F = P = N = that is stable for sufficiently small 
positive values of x '■ if too little energy is provided the laser remains switched off. The 
linear stability of the trivial solution with respect to a small perturbation (F, P, N) is 
given by 



d_ 
dt 



C 1 

x -1 







(9) 



and depends also on the perturbation itself. However, one can see that iV is decoupled 
from F, P and it is always damped. It is convenient to expand both F and P in terms 
of the eigenfunctions A n of C. By projecting the equations for F and P onto the A n , we 
obtain an infinite block-diagonal matrix with blocks of dimension two. The eigenvalues of 
these blocks are given by 



A ra — 



-l + (3 n ± v /(l + /3„) 2 + 4 x 



(10) 
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and the corresponding eigenvectors are 



'i + aA 



1 



V 



X 




A n (x,y) 



(11) 



/ 



where Af n is a normalisation factor chosen so that (v n , u n ) 
sponding eigenvector of the adjoint matrix, 



1. Here v n is the corre- 



1 





(12) 



and the inner product is defined by 



(v,w) = J v T udxdy. 



We have, therefore, that 



A/1 2 



;i + a„) 2 + x- 



(13) 



(14) 



From equation (10) we can calculate that the trivial solution loses its stability through a 
Hopf bifurcation of frequency 



UJf, = 



1 - K(/3 fl ) 



(15) 



when x = Xo(^), where 



Xo (n) = (l + uj 2 n ) 



(16) 



The index n is defined as an element of the set J C N of those values of n that minimise 
Xo( n ) over & H the integers. 

In the presence of unperturbed symmetry there may be degenerate eigenvalues. In other 
words, the set J has more than one element. Eigenvectors corresponding to eigenvalues 
A ra ,n G J are called active modes. In general, we can assume that the number M of 
elements of J is finite, in which case we say that the unperturbed system has degeneracy 
M: all the active modes have eigenvalue with zero real part and imaginary part given 
by (15). All the other modes have eigenvalue with negative real part. If the system is 
slightly perturbed the degeneracy is either partially or totally lifted and the eigenvalues 
of the active modes are no longer all equal. However, we assume that the perturbation is 
sufficiently small and that x is sufficiently close to Xo(^) to ensure that 



£(An) < Vn ^ J. 



(17) 



5 



In other words, the modes that do not belong to the set J are still linearly damped. We 
can separate the field u = (F, P, N) T into components along the space spanned by the 
active modes and its orthogonal complement: 

u = Yl f"< u ™ + u ^ ( 18 ) 

n<=J 

where u ± = (F±, P±, N) T is defined by 

(v„,Ul)=0, Vnej. (19) 

The equations for the dynamics of the active modes can be obtained by a centre man- 
ifold reduction [32] of equations (3-5). The nonlinearities in these latter equations are 
quadratic. As a consequence, the centre manifold reduction is greatly simplified and the 
equations for the active modes can be obtained up to third order by "simply" projecting 
the equations (3-5) onto the active modes: 

d „ / du\ 

7t f " =("»•*)• neJ ' (20) 

After rather lengthy calculations (detailed in Appendix A) we obtain that the amplitudes 
of the active modes are given by 

4/n = KJn ~ X J2 B n-Vjk (A», Aj A k A p ) fj f k f p , 71 € J (21) 

where 

(l + A p )[l + (A, + A fc )/2] 
^npjk n i / \ i \ \ a r \ r \ r iTr ■ \ ) 



[1 + (\j + X^hWnN^Nk ' 



3 Two modes and D 2 geometrical symmetry: an example 



3. 1 Introduction 



Equation (21) is very general: it can be used to describe the bifurcation structure of a mean 
field limit two level Maxwell-Bloch laser at an eigenvalue of arbitrary (finite) multiplicity. 
The methods used in deriving it are also completely general. However, the great advantage 
of using a two level Maxwell-Bloch laser model is that the nonlinear terms are relatively 
simple and hence the centre manifold reduction of the laser equations contains fairly few 
and manageable terms. 

In order to illustrate the use of equation (21) we consider in this section a particularly 
simple case: we assume that the perfectly symmetric laser has two degenerate modes at 
threshold and that symmetry breaking, for example as caused by astigmatism, reduces 
the geometrical symmetry from 0(2) to D 2 and lifts the degeneracy by changing the mode 
resonance frequency. The geometrical symmetry of the astigmatic cavity is combined with 
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the phase invariance S 1 - symmetry of the laser equations to give the symmetry group 
Z 2 x S 1 (see below). This problem has been studied in depth by many authors [14,16,17]. 
Dangelmayr and Knobloch [14] have provided a fairly comprehensive analysis of all the 
possible bifurcation scenarios induced by the symmetry breaking from 0(2) x S* 1 to Z2 x S 1 
in a generic system in normal form. In the laser context, various authors [9,10,29,33-35] 
have studied the effect of symmetry breaking on mode interactions for simple laser modes. 
In a similar vein, Vladimirov [36] studies the interaction between longitudinal modes in 
a bidirectional ring laser. 

This section contributes to these studies by highlighting the fact that the two interacting 
modes need not be standard Gauss-Hermite modes but can be any cavity mode: for 
example they could be the modes of a waveguide cavity as studied in [6-8] . Moreover, we 
include in our analysis a bifurcation scenario with unexpected additional symmetry that 
is excluded from the generic studies in [14] and which arises when the frequency ujq of 
the unperturbed cavity mode and the frequency uja of atomic transmission coincide. We 
now explain this in more detail. 

Perturbations that break the symmetry typically cause the two modes to have slightly 
different frequencies. By a suitable choice of reference frequency we can write their coef- 
ficients j3 n as 

P n = K {-l + i[8 + (-l) n ri]} n = l,2 (23) 

where k is the mode amplitude decay rate and 8 — 8/k is the scaled cavity-atom detuning 
parameter. We have relabelled the modes so that the two active modes have indices 1 and 
2 and we have changed the reference frequency so that the active modes have frequencies 
k(8 ±77). In this notation 77 plays the role of symmetry breaking parameter: if 77 = 
the two modes are degenerate and the system undergoes a standard 0(2) x S* 1 symmetric 
Hopf bifurcation. If 8 — ±77 then atomic transition is resonant with mode n — 1 and n = 2 
respectively. The mode with frequency closest to the atomic frequency (so that S(/3 n ) is 
closest to zero) is the first to become unstable as the pump parameter is increased from 
below to above threshold. However, if the atomic frequency is exactly in between those of 
the two modes, i.e. if 8 = 0, then the laser has no criterion by which to choose the most 
unstable mode and like Buridan's ass does not know what to do. Unlike the unfortunate 
animal, though, it does not die of hunger but exhibits dynamical behaviour that is an 
interesting "combination" of features of the dynamics observed for 8 small and nonzero. 

3.2 Centre Manifold Equations 
We write the pump parameter x as 

X = «(l + e). (24) 

From (23) and (10) we see that at resonance 8 = ±77 the threshold value of the pump 
parameter is x = K i that is e = 0. The assumption of being close to threshold therefore 
becomes that e is a small parameter. Moreover, we assume also |5±?7| <C 1. To write the 
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equations for the two modal amplitudes we simplify (21) by first expanding the coefficient 
of the linear term up to first order in e and up to second order in 5 ± r\ : 



A n = us - fi 3 [5 + (-l) n 7?] 2 + in[5 + (-l) n r)} + ■■■ (25) 

where 

(26) 

This expansion ensures that we keep track of all the physically relevant features of the 
model: in particular, the second order term in S±i] is responsible for making the threshold 
value of the pump parameter dependent on the detuning. We also expand the coefficients 
of the nonlinear terms, but only up to order zero. In this limit equation (21) reduces to 



j t f n =f,{e-^[5+ (-l)^f + i[S + (-l) n v]}fn 

~Y^- K E (A,, A 3 A k A p )fJ k f p , n = 1, 2 . (27) 

j,k,p=i 

In order to determine which of the nonlinear terms of (27) are nonzero we use the fact that 
the geometrical symmetry of the unperturbed system is 0(2) and that the perturbation 
is small so that we can compute the projection integrals (A n , AjA k A p ) assuming that the 
modes A n (x,y) are the modes of the system with full symmetry. 

The (complex) 2-dimensional eigenspace for the linearisation of (27) at the first threshold, 
i.e. at the bifurcation point where the zero field solution loses its stability, must be a 
subspace of some particular isotypic component of the action of 0(2) on the L 2 function 
space 7i spanned by the cavity modes (see [3, Theorem XII, 3.5] for example), and we 
can therefore assume that the modes A n (x,y), represented in polar coordinates (r, iji) as 
A n (r, if)), can be written as (complex) linear combinations of cos(?ra/>) and sm(rmp) for 
some to e N. If we suppose that the x and y axis are aligned with the symmetry axis of 
the perturbed system then we may take 

Ai(r, 1])) — A (r) cos(to/0) and A 2 (r, ip) = A (r) sm(mip), (28) 

where A (r) represents the radial profile of the two modes. The projection integrals in (27) 
have values 

37r/4 n = j = k = p 



(A n , AjA k A p ) = M < 



n=j^k=p 
7r/4 [all permutations of the (29) 

four indices are allowed] 

otherwise 
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where the constant M. is given by 

M= Al(r)Al(r)rdr. 
Jo 

Using (29) we can write (27) as 



(30) 



-f n =f,{e-^[5+ (-l)>f + i[S + (-!)>]}/„ 



(3\f n \ 2 + 2\f p \ 2 )fn + Lfl , n = 1, 2, p = 3 - n. 



4(1 + re) 

To simplify the notation we scale time and the modal amplitudes by defining 



r = fit, a + 
so that we can rewrite (31) as 



4(1 + re) 



/i, 



a. 



ttM 



4(1 + re) 



/ 2 , 



£ — fi 2 (5 — i]) 2 + i(5 — if) a + — ^3|a + | 2 + 2|«_| 2 ) a + — a + a 2 _ , 
e — ji 2 (5 + r/) 2 + i{5 + rj) a_ — (2|a + | 2 + 3|a_| 2 ^) a_ — a 2 + a_ . 



(31) 



(32) 



(33) 
(34) 



These are the expansions to third order of the equations (3-5) reduced to the centre 
manifold, and are the equations that we study for the remainder of this paper. 



3. 3 Symmetries 



The system (33, 34) is unchanged by sign-reversal of either a + or ct_, i.e. it is equivariant 
with respect to the symmetry group D 2 = Z 2 (cr x ) x Z 2 ((Jy), where we use the notation 
Zr(Q to denote the cyclic group of order r generated by the element £. The actions of the 
group generators are: 



(35) 



and they represent the geometrical symmetries of the laser. Using (28) we see that, for 
example, a x corresponds to a reflection with respect to the x-axis: this leaves the cosine 
mode (amplitude a + ) unaltered and changes the sign of the sine mode (amplitude 

The system (33, 34) is also equivariant with respect to the S 1 action corresponding to 
global change of the phases of the two amplitudes: 











\ a -J 













<j y : 

















R 



v ■ 



UA 












\ a -j 







(36) 
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In other words, the amplitudes a + and a_ are defined only up to a common phase factor. 
Therefore the total symmetry of equations (33, 34) is the combination of the geometrical 
and the phase invariance symmetry. However, we note that we can write 

a x — o y o R n , (37) 

so the total symmetry group can be taken to be Z 2 x S 1 , where Z 2 is generated by one of 
either o x or o y . The other geometrical symmetry can be obtained using (37). Therefore, 
the geometrical symmetry group D 2 is recovered as a subgroup of Z 2 x S 1 . 

We are now in a position to describe the bifurcation structure and interaction of the two 
modes following the extensive analysis in [14]. We summarise here the main results as 
they apply to our case. In the notation of [14] we have a = b = — 1 which places us on the 
lower boundary of region III(c) in Figure 4 of [14], from which the theoretical bifurcation 
diagrams can then be deduced using Figures 6 and 7 of [14]. In our Figure 1 we give 
versions of these diagrams as observed numerically in our context. The next subsection is 
devoted to explaining this Figure. 



34 The case 5^0. 



We consider first the case 5 > 0; the results obtained are also valid for 5 < after 
exchanging a + with a_. The important special case 8 = of cavity resonance, excluded 
from consideration in [14], will be discussed in the next subsection. 

The phase invariance of (33) and (34) implied by the S 1 symmetry makes the bifurcation 
analysis rather cumbersome. Therefore, it is convenient as in [14] to introduce the variables 
A, <3> and 6 defined by: 



v = Acos($/2)e iei , 
w= Asin($/2)e i02 , 

e= 0i + 2 , 



(38) 
(39) 
(40) 



where 



v = a + — ia- and w = a + — ia- 



(41) 



We can reduce the two complex equations (33) and (34) to three real equations for A, $ 
and : 



dA 
~dl 

dQ 
~dt 



e-/j 2 (5 2 + rf 



A - 



1 9 

1 + - sin $ 
2 



A 3 + cA sin $ cos cos a 



— — sin(2$) + 2c [cos cos a cos $ — sin sin a] 

Q 

—2 [sin cos a + cos sin a cos $] 

sin $ 



(42) 
(43) 

(44) 



10 



[cf. [14]) where c is defined by 



ce ia = r](2fi 2 5 - i) 



(45) 



and plays the role of a (complex) symmetry breaking parameter in that it is zero if and 
only if the 0(2) symmetry of the system is present. Note that in terms of A, $ and 6 
both a x and a y have the same realisation 



(46) 



thus highlighting the fact that the symmetry of the problem can ultimately be considered 
to be Z 2 , once the phase invariance S 1 symmetry is factored out. 

At e — E\ = fi 2 (5 — rf) 2 the trivial solution a± = (A = = 0, $ = 7r/2) loses its 
stability through a Hopf bifurcation (bifurcation point Pi in Figure 1). A new nonzero 
single mode solution (standing wave solution SW in the notation of [14] and in Figure 1) 
with frequency 





(A 




( A \ 








IT - $ 








{ - e ) 



w+ = 6-ri, 



and amplitudes 



I i o 1 / 2 2 \ 

l«+l = o{ e ~ ^ u +) ' 



a_ = 



(47) 



(48) 



appears and is stable. Note that in the context of equations (42-44) this bifurcation is a 
pitchfork and the new stable solution is stationary with values 



' 2 2 



) , $ = tt/2, = 0. 



(49) 



This simplification is a result of having factored out the S 1 (phase invariance) group. The 
stationary point is a focus for 



e < e N = 6i] + fi 2 (5 - i]f 



(50) 



and is a node otherwise. 



There is an additional bifurcation from the zero solution. At e = e e = fi 2 (5 + rf) 2 there is 
a Hopf bifurcation to the unstable solution (SW n in the notation of [14]) 



a + = 0, 



\a. 



e- n 2 (5 + r)f 



(51) 



In terms of A, $ and this is a pitchfork bifurcation to the unstable solution 



A 2 = - 



e-fi 2 (5 + r]) 2 , $ = tt/2, = tt. 



(52) 



11 




4 
3 
2 
1 




5 6 



A 




tw ; 


_fii 


^2 


e' 



2 4 6 8 10 




2 4 6 8 10 




4 
3 
2 
1 




A 





















2 4 


6 8 10 



Fig. 1. Bifurcation diagrams for k = 1 and rj = 1 (which imply that 82 = 2) and different 
values of the detuning parameter: 5 = {1,1.9,2.1,3} from (A) to (D). The solid (dashed) lines 
correspond to stable (unstable) solutions. We have not drawn the unstable branch SW n . The 
bifurcation points Pj, j = 1, ..,5 have horizontal coordinate £j. The vertical coordinate of the 
branches of stationary solutions is their amplitude A. That of the Hopf branch is the maximum 
of the amplitude during one period. For this reason in panel (A) the Hopf branch does not touch 
the limit point. 

There are no further bifurcations on this branch. 

The SW branch may have a secondary Hopf bifurcation to a two-frequency solution (H 
in Figure 1) as the bifurcation parameter is increased to the value 



e 2 = fi 



(5 + ri) 2 + 85r) 



(53) 



This solution is characterised by both modes having non-zero amplitude, but different 
frequency so that beating between the two modes occurs. 

As 5 tends to zero 62 tends to e±, so that the secondary Hopf instability occurs closer and 
closer to the primary bifurcation from the trivial solution (see panel (A) in Figure 1). 
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Fig. 2. (A) Blue sky bifurcation that ends the existence of the periodic orbit in panel (A) of 
Figure 1 at e = 5.399. Two saddle node bifurcations (circle: stable fixed point, square: unstable 
fixed point) appear at symmetrical locations along the periodic orbit. (B) Disappearance of the 
periodic orbit in the case of panel (B) of Figure 1 at e = 6.098. The periodic orbit clashes against 
the heteroclinic connections between the two unstable fixed points. All the parameters are as in 
Figure 1. 



This secondary Hopf bifurcation point P 2 exists only in the range 

1 

2^ : 



< S < 5 2 = — . (54) 



In fact, another two stationary solutions appear at a third bifurcation point P 3 on the 
SW branch, at 



3r/ 



e = e 3 = ^ + n \q z + 5 2 + A5r]\ . (55) 

In these solutions the two modes are frequency locked to each other and thus the intensity 
pattern is stationary. 

If 5 — 82 then e 2 = £3 and the two bifurcation points P 2 and P3 collide. The solutions 
that appear at P3 (TW' in Figure 1) have a + and a_ both nonzero ($ 7^ 7r/2 and 0^0) 
and are interchanged by the symmetry operation a defined by (46). The bifurcation is 
sub-critical if 




0<5<5 4 = \ -6 2 , (56) 



(see panels (B) and (C) of Figure 1) and has a limit point P4 at 

£ = £4 = 2v / 6~7/ + /iV + 5 2 ). (57) 

If 5 = 64 then e 3 = e 4 and the limit point P4 hits the SW branch at P3. For larger values 
of S the bifurcation is supercritical (see panel (D) of Figure 1). 

The secondary Hopf branch ends either at P4 or, for larger values of 5, at a point P 5 
on the the unstable branch that joins P3 and P4 (see panels (B) and (C) of Figure 1). 
In the former case the periodic orbit disappears in a blue sky bifurcation (panel (A) of 
Figure 2), that is the period tends to infinity until a saddle-node bifurcation creates a 
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Fig. 3. State diagram of the equations (42~44) f° r n = 1 and rj = 1. The line marked Pj 
represents the bifurcation point Pj of Figure 1 as a function of the parameters e and 5. The 
dashed line represents the curve e^: below it the SW point is a focus, above it is a node. The 
large dot next to the letters TB indicates the Takens-Bogdanov bifurcation point. The top part 
of the figure is an enlargement of the rectangular dashed section in the bottom part. The orbits 
drawn are just indicative of the type of dynamics expected in the corresponding parameter region 
and are not intended to be in any particular plane of the A, and space. 

sink and a saddle that 'block' the periodic orbit: compare [37, Fig.3]. In the latter case 
the period also tends to infinity, but the orbit collides with the heteroclinic orbits that 
join the unstable fixed points in the TW branch (see panel (B) of Figure 2). There is no 
analytical expression for the coordinates of P 5 and we have used AUTO [38,39] to plot it 
in Figure 1. Finally, note that even though the point P3 tends to infinity as 5 tends to 
zero, P4 is always finite: therefore the periodic orbit exists for only a finite range of values 
of e. These results are summarised in the state diagram in Figure 3. 
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The state diagram is organised around the bifurcation point of highest codimension (two), 
namely the Z 2 -symmetric Takens-Bogdanov point (indicated by TB in Figure 3) which 
occurs when the points P2, P3 and P5 (see panels A and B of Figure 1) coincide, i.e. at 
5 = 62- In [14] the normal forms in the neighbourhood of this bifurcation point are derived 

as 



in suitable coordinates £ and 77. The symbol 0(5) represents terms of order at least five in £ 
and 77 and the coefficients Mtb and Ntb depend on the parameters of equations (33, 34). 
The condition for a (non degenerate) Takens-Bogdanov singularity is MtbNtb 7^ 0. 
Using the calculations given in Appendix A of [14] we find in our case M TB = 1/3 and 
Ntb = =f25/(18a/2), where the signs refer respectively to the cases 5 > and 5 < 0. 
This implies that the bifurcation geometry in a neighbourhood of the Takens-Bogdanov 
point in Figure 3 is topologically equivalent to that obtained from a versal unfolding of a 
standard model as described for example in [18]. 

The bifurcation geometry in the neighbourhood of a ^-symmetric Takens-Bogdanov point 
is well documented: initial studies by Takens [40] with further analysis given by Carr [32], 
Guckenheimer and Holmes [41], Arnol'd [42], Khorozov [43] (see also [44]) and others show 
that a versal unfolding of the vector field at this singularity gives a bifurcation diagram 
determined by a smooth curve passing through the Takens-Bogdanov point and two rays 
emanating from it not tangent to the curve or each other (in contrast to the non-symmetric 
case). This is indeed the description of the bifurcation diagram in the neighbourhood of 
the point TB in Figure 3. 

In applications, a Takens-Bogdanov bifurcation, whether symmetric or not, is often ac- 
companied by additional characteristic features that arise from interaction with further 
equilibrium states. This can be seen, for example, in the forced van der Pol equation [45,46] 
(see also [41]) and in the dynamics of bulk liquid crystals in a shear flow [47]. The expla- 
nation for these accompanying patterns lies in the geometry of unfolding of a singularity 
of higher degeneracy and thus higher codimension: although it does not appear among the 
vector fields studied, its presence behind the scenes organises the configuration of lower 
codimension bifurcations. In the case of the van der Pol oscillator there is a triangular 
region in the neighbourhood of the (non symmetric) Takens-Bogdanov point (see Figure 
2.1.2 of [41]) whose footprint can recognised as part of the codimension-3 non symmetric 
bifurcations studied by Dumortier et al. [48]. 

An analogous phenomenon occurs here. The organising centre is a degenerate Z 2 -symmetric 
Takens-Bogdanov singularity involving coalescence of five equilibrium states (instead of 
three for the non-symmetric case) and arising when Mtb — in equation (59). As in the 
van der Pol case, there is a 'triangular' region (shaded area in Figure 3) whose features are 
determined by the interactions of the five fixed points. However, we are not aware of a rig- 
orous treatment of this case even though many authors have touched on aspects of it. For 
example, interactions of the five equilibria lead to interesting dynamical phenomena in cer- 
tain fluid dynamical convection problems [12]. In the case of lasers, Lopez- Ruiz et al. [33] 



£ = r/ + 0(5), 

7) = M TB £ 3 + iV Ti? £ 2 77 + 0(5), 



(58) 
(59) 
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Fig. 4. Representation of the periodic orbits on the equations (42-44) i n (Q,^)- spherical coor- 
dinates for different values of the detuning parameter 5. The orbits are contained in the nearer 
(further) hemisphere if 5 > (5 < 0). The limiting case 5 = corresponds to an orbit co- 
inciding with a meridian of the sphere (thick line). In these simulations k = 1, r\ = 1 and 
5 = {0.75, 0.50, 0.25, 0.00, -0.25, -0.50, -0.75}. 

give a bifurcation analysis only in the neighbourhood of the symmetric Takens-Bogdanov 
point. Krauskopf et al. [37,49,11,35] give spectacular and informative bifurcation diagrams 
away from the Takens-Bogdanov point, for a model that does not, however, involve the 
^-symmetry. Rucklidge [13] considers bifurcations only near the Takens-Bogdanov point 
and with the sign of M TB opposite to ours. 



3.5 The case 5 = 
3. 5. 1 Introduction 

If 5 = the symmetry-breaking parameter c defined by equation (45) is purely imaginary 
and equal to — irj. The physical interpretation is that from the point of view of the laser 
the modes are distinct because they have different resonant frequencies (77 7^ 0), but also 
there is no reason for the laser to chose one rather than the other as the pump parameter 
is increased through threshold. If instead 5^0 then one mode has a lower threshold than 
the other, reflected in the fact that the symmetry-breaking parameter has a non-zero real 
part. 

The solutions of equations (42-44) can be conveniently represented in spherical coordi- 
nates: the variable G plays the role of the polar angle and the variable $ that of the 
azimuthal angle (see Figure 4). The two fixed points A 7^ 0, $ = n/2 and 6 = {0,7r} 
correspond to diametrically opposed points in the equatorial plane and the periodic orbits 
are closed loops around them. Moreover, the sign of the detuning 5 determines in which 
hemisphere the periodic orbit is situated. The orbits for opposite values of 5 are related 
by the symmetry a defined by (46). 

If 5 = the orbit is constrained to the meridian O = 71/ 2 that separates the two hemi- 
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spheres (see Figure 4) and thus has higher symmetry. This suggests in turn that the 
model may have an additional symmetry if 5 = 0. Therefore if we identify this additional 
symmetry we can use the tools of symmetric bifurcation theory [3,4] to analyse the 5 = 
case and its response to symmetry-breaking perturbations. This is the approach we follow 
throughout the rest of this section. 



3.5.2 Additional discrete symmetry 

As a first step we identify the symmetry group of the equations (33, 34) when 5 = 0. In 
this case the equations (33, 34) become 



a + = (A — irj)a + — (3|a + | 2 + 2|a_| 2 )a + — a + ai 
a_ = (A + ir))ot- — (2|a + | 2 + 3|a_| 2 )o;_ — a^a 2 



(60) 



where A = e — fi 2 i] 2 . Under the coordinate change (41) these correspond to the perturbed 
equations (2.7) of [14] 



v = g(X, \v\ 2 , \w\ 2 )v + cw 
w = g(X, \v\ 2 , \w\ 2 )w + cv 



(61) 



which have Z 2 x S 1 symmetry given by 





and 




-up 



(62) 



w 



as shown earlier using equation (37). However, if 5 = the perturbation coefficient c is 
purely imaginary, and this case is excluded from consideration in [14]. Since g is real there 
is now an additional symmetry of order 2: 



or, equivalently, 



r : 



t : 



(')' 


-> —i 
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(*-\ 

i — > 




\a + ) 



(63) 



(64) 



Combined with the S 1 action 









\ a -J 




\ a -l 



(65) 
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• Re(a + ) 



Fig. 5. Schematic representation of the intensity patterns in Fix(Si), equation (71). The black 
areas represent regions of high intensity. Their position along the thin circle indicates the cor- 
responding value of the argument of a + . 

this generates an action of 0(2). Moreover, the symmetry 



P ■ 













ia + 








\ a -J 




y —ia^ j 



(66) 



of order 4 commutes with both r and R v . This therefore gives as total symmetry group 
the subgroup Z 4 x 0(2) of the group S l x 0(2) generated by 









\ a -J 







(67) 



and 0(2) as above. Notice that the S 1 x 0(2) action here corresponds to that in [14] under 
direct substitution of v,w in place of a + ,o;_ although in our problem these coordinates 
are in fact related by (41). We shall return to this significant point later. 

Observe that there is redundancy in the group action since the element ( = p 2 R n acts 
trivially. Removing this redundancy yields the quotient group (Z± x 0(2))/Z 2 (C) which 
is cumbersome to work with both conceptually and notationally - in contrast to the 
(D 2 x S l )/Z 2 (a x ) case handled earlier which reduced simply to Z 2 x S 1 . We therefore 
prefer to focus on the group Z4 x 0(2) with its transparent algebraic structure, while 
recognising that its action on C 2 = R 4 is not faithful. 



Once the symmetry group of the dynamical system has been established, we proceed to 
identify its subgroups and their fixed point spaces. The fixed point space Fix(E) of a given 
subgroup £ is the vector space whose points are left unchanged by the action of E . The 
importance of fixed point spaces for bifurcation analysis is that they are invariant under 
the dynamics. In particular if Fix(E) is one-dimensional then typically at bifurcation there 
must be a branch of stationary solutions lying in Fix(E) (Equivariant Branching Lemma) 
and thus exhibiting the symmetries E . If Fix(E) is two-dimensional and the eigenvalues 
at bifurcation are purely imaginary then typically there is a branch of periodic solutions 
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lying in that subspace (Equivariant Hopf Theorem). See [3] or [4] for precise statements 
of these results. 



The subgroups of Z 4 x 0(2) with 2-dimensional fixed point subspaces are found among 
those for S* 1 x 0(2) (see [3]) and are (up to conjugacy) as follows: 



Ei : Z 2 {t) x Z 2 (p 2 R 7T ), E 2 : Z A (pR-«). 
Their respective fixed point spaces are 



(68) 



Fix(Sx) 




Fix(E 2 ) 




(69) 



We next identify these fixed point spaces in terms of laser patterns. We know from equa- 
tion (28) that a + and a_ are the amplitudes of a cosine and sine mode respectively. In 
particular, for the purpose of this example we can assume that m — 1 in (28) so that the 
expansion of the electric field can be written as 



F(r, tp, t) = a + (t)A Q (r) cos(^) + a-(t)A (r) sin(^). 
The intensity pattern that corresponds to the elements of Fix(Ei) is 

I^FixCSi)! 2 = |A (r)|>+| 2 {1 + cos [2arg(a+)] sin(2^)} . 



(70) 



(71) 



As the argument of a + goes from to 2tt the intensity pattern changes according to 
Figure 5. The elements of Fix(E2) have an intensity pattern given by 



|^Fix(E 2 ) | 2 = |A (r)|> + | 2 cos 2 (^). 



(72) 



This corresponds to two bright spots aligned along the horizontal axis, independently of 
the phase of a + . There is also the corresponding sine pattern: this corresponds to Fix(E' 2 ) 
where E 2 is the conjugate E 2 = tE 2 t~ 1 of E 2 . 



Given a stationary or periodic solution of a symmetric dynamical system we obtain other 
solutions by applying the elements of the symmetry group to it. It is important therefore, 
to understand how the full group Z 4 x 0(2) acts on the fixed point spaces Ei and E 2 . 

Let a denote the vector (a+, aJ) T . If ^ a e Fix(Ei) then pa e Fix(Ei) as p commutes 
with r. On the other hand, R^a e Fix(Ei) just when tp — 0, it. Therefore the Z4 x 0(2) 
orbit of a consists of two circles intersecting Fix(Ex) at {a, — a} and {pa, —pa}. In 
physical terms this result means that if there is a solution in Fix(Ex), then there are 
infinitely many copies of it, each obtained by multiplication by a phase factor. However, 
if the phase factor is equal to 7r then the net effect is not to obtain a different solution 
but to have moved half a period in Fix(Ex). 

If ^ a e Fix(E 2 ) then pa e Fix(E 2 ) and R^a e Fix(E 2 ) so the Z 4 x 0(2) orbit of 
a consists of two circles C and tC lying in Fix(E 2 ) and Fix(E' 2 ) respectively. In other 
words, the net effect of the group action is to multiply a sine or cosine solution by a phase 
factor. 



19 



We can combine these results and make use of the Equivariant Hopf Theorem [3,4] to 
deduce the existence of branches of periodic solutions with these symmetries that bifurcate 
from the trivial solution. 

Theorem 1 For the equations (60) there exist (at least) the following periodic solutions 
branching from the trivial solution as A passes through zero: 



(1) a family of solutions of the form 



\r x e 



i<P e —i(Tft+0) 



(73) 



where r\ ~ A 2 and cp and 9 are arbitrary. Here W v has reflection symmetry under 
R^tR' 1 = R 2lfi r and a Z 4 spatio-temporal symmetry under which increasing time by 
7r/ (217) is equivalent to applying p ; 
(2) a pair of solutions of the form 



V + ,V_ 



( r - xe i{vt+e) 



\ 







K r x e 





-i(r)t+0) 



(74) 



respectively, where r\ ~ A 2 and 9 is arbitrary. Each solution has spatio-temporal 
symmetry SO (2) given by shifting time by At and rotating by 9 = —r]At. 

Under time evolution the solution W v cycles across all the patterns in Figure 5. The 
solutions V + and V- correspond to the solutions SW and SW n , equations (48) and (51). 
We will see below that for our particular system (60) these are the only solutions that 
branch from the trivial solution and, moreover, that the W v solution is stable and the V± 
are unstable. 



3.5.3 Disguised rotational symmetry 

Interestingly as it turns out, the relations between the solutions for 5 = and those for 
5^0 described above can be elucidated by rewriting equations (60) in terms of new and 
somewhat artificial variables that restore to the equations a higher degree of symmetry. 
The linear terms of (60) commute with rotational symmetries 



(75) 



e^ 2 a. 



for arbitrary </?i, </? 2 , and so Takens' theory of normal forms [18] suggests the possibility of 
using a nonlinear coordinate change to remove (to any desired order) all nonlinear terms 
that do not share this symmetry. It is important to stress that the change of variables is 
nonlinear: although the linear transformation (41) removes the non-0 (2)-symmetric cubic 
terms a + a 2 _ and a 2 + a^ from equations (60), the linear terms in v and w in the resulting 
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equations (61) at the same time lose their full rotational symmetry. Normal form theory 
states that the nonlinear transformation will succeed when there are no resonances (that 
is, linear relations with integer coefficients) between the eigenvalues at the bifurcation 
point. In our case the eigenvalues at the bifurcation point are rkir] and there are of course 
strong resonances, but careful inspection of the normal form method shows that the terms 
a + a 2 _ and a_a\ can nevertheless be removed. To be precise: 

Theorem 2 The nonlinear coordinate transformation 









I- 1 


f M 2 -\ 


\ a -J 






1 2A + Air] 





(76) 



converts the equations (60) into 

$ + = (A - irj)(3 + - (3|/3 + | 2 + 2\(3.\ 2 )(3 + 
/3_ = (A + ^_-(2|/3 + | 2 + 3|/3_| 2 )/5_ 

up to order 3. 

Proof. Direct verification, or appeal to general theory [18] and the observation that it is 
the difference and not the sum of the eigenvalues A ± irj whose vanishing provides the 
obstruction to removing the non-rotational cubic terms. 

Observe that the equations (77) have the form (61) with e = and with v, w corresponding 
to (3+,f3- , in contrast to the original correspondence defined by (76) and (41). Using the 
framework of [14] we can therefore express the symmetry as follows: 

Corollary 1 After changing to (3 = (f3 + ,/3-) T coordinates and truncating at order 3 the 
system (60) has S 1 x 0(2) symmetry with 0(2) generated by r, R v and with S l -action Pg. 

Remark 1. We recover the independent phase symmetries (75) by taking ipx — ip + 9 and 

ip 2 = ip - 6 . 

Remark 2. If in (60) the two complex eigenvalues had been A + in (twice) rather than 
A ± if] then the equations would have D A x S* 1 symmetry and (60) would already be in 
normal form: the terms (a + a 2 _ , a^a 2 + ) could not be removed. Compare [50], Sect. 3.1. 

Remark 3. We have been working with Taylor series truncations at order 3. In fact, normal 
form theory would allow us to transform away all those higher order terms in a full Taylor 
series extension of (60) that did not have S* 1 x 0(2) symmetry. However, the higher the 
order the smaller the likely region of validity, and there is no guarantee that the system 
itself exhibits strict S* 1 x 0(2) symmetry. This does not affect robust (structurally stable) 
features of the bifurcation scenario, but can affect sensitive features such as heteroclinic 
connections and families of orbits on an invariant torus. See [50] for a careful discussion 
of this issue for a system closely related to ours. 
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Fig. 6. Changes in the primary bifurcation point as 5 is increased from zero. The letters outside 
the brackets identify the branches with those in Figure 1 and refer to the equations in terms 
of a± . The letters in square brackets give the names of the branches as viewed in terms of the 
variables (v, w). 

The advantage of using the variables (3 is that we can now recast the problem of connecting 
the solutions for 5 — to those for 5 ^ as the problem of breaking the 0(2) symmetry 
studied in [14]. There it is shown that the fully symmetric system has a standing wave 
solution v — w and two travelling wave solutions (v , 0) and (0,w). We indicate these 
standing and travelling wave solutions with the symbols SW* and TW* respectively. 
Translating in terms of a + and a_ by identifying (v,w) with (f3 + , f3_) and noting that 



we thus recover the results of Theorem 1 together with the additional information: 

Corollary 2 The transformation oc i— > (3 converts the periodic solutions with Z 4 spatio- 
temporal symmetry into rotationally symmetric ones, at least up to third order. 

In particular the Hopf solution W v (denoted by H in Figure 1) corresponds to SW*, while 
the two solutions V± (SW in Figure 1) correspond to the two solutions TW*. 

Furthermore, it is shown in [14] that as the real part of the symmetry breaking parameter 
e in equations (61) moves away from zero the two TW* branches separate and the SW* 
solution appears as a secondary bifurcation from one of the TW* branches. In terms of 
the variables a this implies that as 5 moves away from zero the two SW branches have 
different starting points, one at e — E\ the other at e = e§. Moreover the Hopf branch 
becomes a secondary branch off the first SW branch. This is illustrated in Figure 6. 

This picture is confirmed by standard bifurcation analysis for systems with two purely 
imaginary eigenvalues [18,41]. Using the (3 coordinates for both S zero and nonzero (but 
small) the parameters e, S provide versal unfolding parameters about the organising centre 
(e, 5) = (n 2 rf, 0). In particular, we can see from Figure 7.5.2 of [41] (redrawn here in terms 
of e and 8 as Figure 7) how the branching of the solutions W v and V± that is simultaneous 
when 5 = changes when 5 ^ 0: a succession of two primary bifurcations with either 
a + = or a- = is followed by a secondary (Hopf) bifurcation to a torus of periodic 
solutions. From (73) we also see that the period of these solutions is approximately given 
by 2n/r]. 



a + = ^ (3 + = 0, a- = (3. 



0, a + = a- 4=> (3 + — (3- , 



(78) 
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Fig. 7. Bifurcation set in the neighbourhood of (e,5) = (fj, 2 w 2 ,0) in terms of the variables 0± 
(adapted from [41 J, Figure 7.5.2 case la). The axes of the small inserts are \(3±\ and the filled 
and open circles indicate stable and unstable fixed points respectively. 

To conclude the analysis of the bifurcation diagram for 5 = we note that the limit points 
of the TW branches of panel A of Figure 1 exist for all sufficiently small \8\ (including 
5 = 0). Exactly as in the case 5^0 they are created at a symmetric pair of saddle- node 
bifurcations of periodic orbits that annihilate the torus (blue sky). These branches, found 
by direct calculation as in [14], are not connected to either of the a + = , a_ = branches 
when 5 = and so are not detectable by local bifurcation analysis in that case. 



3.6 Numerical verification 



In order to illustrate the range of applicability of the normal forms derived in this article we 
have integrated numerically the equations (3-5) that correspond to the cavity in Figure 8. 
The empty cavity modes can be expressed in terms of Hermite polynomials and are given 
by [51,52]: 



A pq (x, y, z) 



21-P-Q 



^ 7TW x (z)w y (z) p\ q\ p _ 



' V2x ' 






e -x 2 /w 2 x {z) e -y 2 /w 2 y (z) 


w x (z)_ 




_Wy{z)_ 





(79) 



The indices p and q are nonnegative integers: p is the number of zeros of the field along 
the x-axis whilst q is the number of zeros along the y-axis. The functions w x (z) and 
Wy(z) are z dependent scaling factors that represent the field beam waist, i.e. the spot 
size radius, along the x and y axis. They are functions respectively of R x {0) and R y {0) 
and are therefore different for all 9^0. For each value of z the cavity modes form an 
orthonormal basis for the space 7i. 



We have assumed that there are only two active modes, (p, q) = {(1, 0), (0, 1)}. The first 
mode is a "cosine" mode, i.e. its dependence on the polar angle tp is cos(</?), the second 
is a "sine" mode. In order to associate the solutions discussed in the previous section to 
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Fig. 8. Schematic diagram of a ring cavity laser. The z coordinate is along the cavity axis of 
length L and the points z = and z = L coincide. 



laser patterns one should keep in mind that the intensity profiles of the cosine and sine 
mode consists of two bright spots aligned horizontally and vertically respectively. 

The hypothesis that only two modes are active may not be realistic if the laser is operated 
well above threshold, but, on the other hand, the purpose of these simulations is to 
illustrate the application of the normal forms: these are only valid relatively close to 
threshold as they have been derived under the assumption that the field amplitudes are 
small. 

We have scaled the longitudinal coordinate z and the radius of curvature of the mirror 
R to the cavity length L. The transverse dimensions x and y and the beam waists w x (z) 



and Wy(z) are scaled to y LX/n, where A is the laser wave length. In these units the beam 
waists in the plane z = are given by [52]: 



w x (o) 



Rcos6 1 
2 4 



1/4 



Wy(0) 



R 



2 co^e 



1/4 



(80) 



The phase shift per cavity round-trip of the mode (p, q) is given by 

1 



Id. 



p,q 



p + - J cu x + { q+ - j u y , 



(81) 



where 



RcosO 



, uj y = arccos 



cos 9 



(82) 



so that the symmetry breaking parameter i] used in (23) is given by 



V = 



\W X ~ Uy\ 

2(1 -ft)' 



(83) 
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where 1Z represents the total field reflectance of the cavity mirrors. The field decay rate 
k is related to the mirror reflectance and the cavity round trip time T c by 



k = 



l-K 

T 



(84) 



Finally, the scaling factor M. defined in (30) can be obtained by writing in polar coor- 
dinates in the plane z = the cavity modes A pq given by equation (79) assuming that 
w x = w y (i.e. that the symmetry breaking is small enough for its effects to be neglected 
in the coefficients of the nonlinear terms). If we separate the radial and angular part and 
integrate them separately we obtain that 



We have chosen the values of the laser parameters in order to maximise the speed of 
integration: 7 = 1, 1Z — 0.97, T c = 0.03 (hence k = 1). Finally, the radius of curvature of 
the mirror was set to R = 3/2 and the angle 6 to 7r/32 so that the symmetry breaking 
parameter is relatively small, 77 = 0.11377. We have simulated the case of one mode being 
in resonance with the atomic medium by selecting a value of the detuning equal to the 
mode frequency, that is 5 = 77 . 

The program to integrate equations (3-5) represents the polarisation and population inver- 
sion on a rectangular grid thus transforming the partial differential equation into ordinary 
differential equations for the values of the fields at the grid points. These are integrated 
using a variable order variable step method [53]. The projection integral in equation (3) 
is computed using Gaussian quadrature [54]. 

We have run a set of simulations for different values of the pump parameter \ = K (l +£)■ 
The results of these simulations are shown in Figure 9 overlaid on the bifurcation diagram 
of equations (33, 34). For small values of e the agreement between the full equations and 
the predictions of the normal form is excellent. The laser switches on in a single mode 
stationary solution (see inset SW of Figure 9). This becomes unstable for larger pump 
values and a periodic orbit appears that oscillates between two patterns (see inset H of 
Figure 9). As the smallness parameter e becomes larger the agreement between the normal 
forms and the full equations becomes worse. In both cases the periodic orbit disappears 
and is replaced by a stationary two mode solution (panel TW in Figure 9), but the value 
of the pump parameter at which this happens is much larger for the full equations than 
for the normal forms. However, it is quite unlikely that only two modes of the laser would 
be active for these values of the pump parameter and is therefore debatable whether 
one should rely on the normal forms as a guide to the laser behaviour in this region of 
parameter space. Be that as it may, it is clear from the numerical simulations that the 
normal forms are an excellent vehicle for understanding the dynamics of a laser with 
broken symmetry both on a qualitative and quantitative level for values of the smallness 
parameter up to roughly 1/2. 
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Fig. 9. Bifurcation diagram of equations (33, 34) for k = 1 and rj = 0.11377, S = rj. The circles, 
diamonds and squares represent the results of numerical integration of the equations (3-5) for the 
cavity in Figure 8 using the parameters detailed in the text. The circles correspond to stationary 
single mode solutions (SW), the diamonds to periodic solutions (H) and the square represents 
a stationary two mode solution (TW). The images in the inserts show the stationary field in 
the case of the SW and TW' solution and the two oscillating patterns in the case of the periodic 
solution. 

4 Conclusions 

Normal forms have long been known as an extremely useful tool for studying the bifurca- 
tion structure of generic models derived solely on the basis of symmetry considerations. 
However, the connection between the parameters of a specific model and the coefficients 
of the corresponding normal forms is seldom made explicit in the bifurcation literature. In 
this paper we have addressed this issue for the case of a two level Maxwell-Bloch laser and 
have derived the formulae that express the coefficients of the normal forms in terms of the 
parameters of the laser. This is important for two reasons. On the one hand, we can iden- 
tify the regions where the predictions derived from the normal forms are in quantitative 
agreement with the predictions of the full model in contrast to those where the agree- 
ment is only qualitative, as illustrated in Figure 9. On the other hand, we can clearly 
disentangle the "model-independent" features of a theory from the "model-dependent" 
features. Mo del- independent features are, for instance, the type of solutions and the bi- 
furcations observed, which depend on the symmetry and on the first bifurcation and are, 
therefore, common to all models having the same symmetry and the same first bifurca- 
tion. Mo del- dependent features include the positions of the bifurcations and the widths of 
the stability windows in the control parameter space, and depend on important physical 
aspects of each individual model that do not affect the symmetry. In optics, for instance, 
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the number of energy levels included in the description of media interacting with light 
does not change the symmetry of the model. When more than one model is available, the 
ability to separate model-dependent and model-independent features is vital for ascertain- 
ing which model is the most accurate. In principle, this could be done by systematically 
comparing experiments and models over a wide range of control parameters. In practice, 
however, such checks are very difficult to carry out without a knowledge of where to look 
in parameter space for model-dependent features, since outside these regions all models 
with the same symmetry behave in much the same way. 

In addition we have combined in this paper the analysis of two key bifurcation phenomena 
that apply to the Maxwell-Bloch laser and exhibited the relation between them in param- 
eter space. Each of these phenomena has been well-studied in a variety of theoretical and 
applied contexts in the literature, but we are not aware of other studies that incorporate 
the full picture as we have presented it here. The curious interplay of travelling waves and 
standing waves is a new observation. 

The role played by the symmetry generated by r at 8 = in organizing the bifurcations, 
as shown in Figure 6 and 7, is an important new result of section 3. From a physical point 
of view, this additional symmetry is present when the maximum of the gain profile uja is 
between two cavity modes and these modes have the same gain. Up to the second order 
expansion of the gain line in small detunings, this happens for all active media. For large 
detunings, the gain profile needs to be symmetrical with respect to uja, a condition verified 
in most lasers. We can conclude then that this symmetry is useful in understanding the 
interaction of two laser modes for a class of laser much larger than that described by the 
Maxwell-Bloch model used in this paper. 

A final point worth remarking is that our derivation of the normal forms takes place in a 
broad context of breaking rotational symmetry. In this paper we have analysed in detail 
the case of an astigmatic laser, but the same normal forms could be used to study the 
symmetry breaking induced for example by the interaction of a metallic wave guide and 
an aperture, the case studied in [6,7]. 
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A Detailed calculations of the centre manifold expansion 

At the heart of the centre manifold reduction [32] lies the splitting (18) of the fields 
u . Let denote the complex linear subspace of the function space Ti generated by 
the active modes u n , n e J. Centre manifold theory ensures that the local dynamics 
of equations (3-5) can be obtained by considering the restriction of the laser equations 
to the centre manifold of the trivial solution at bifurcation. On this manifold we can 
consider the non-active modes as functions of the active modes, that is u± = uj_(it||). 
The expression for the centre manifold is obtained by substituting this functional relation 
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into the laser equations (3-5). However, this is quite intractable unless we make some 
further simplifications. We therefore assume that we are sufficiently close to threshold 
so that the amplitudes of the active modes are small and we can approximate the centre 
manifold by the leading terms of its power series expansion in the amplitudes of the active 
modes. 

In order to understand which terms of this series we need to compute it is convenient to 
write out in full the equation (20) for the amplitudes of the active modes: 



^ = ^( 1 + A « 10 ) X 

' ~ E f$^-{K, A p )f p - (A n , F ± ) + xJ2 + (A n , P ± ) ^ 



(A.l) 



P ej K 



P ej K 



{A n , A p ) 



P eJ K 



f P -(A n ,P±)+ 



£ ^4r^«> & + WArXp + (A n , (x + N)F ± ) 



"7 



N — -(FP + FP) 



where the symbol {f,g) represents the L2 inner product of two functions f(x,y) and 
g(x,y) . Using the orthogonality relation (19), we can eliminate the terms linear in F± 
and Pj_ in equation (A.l) and write it as 



J t fn = [X (1 + 2A n ) - P n (1 + A n ) 2 ] fn 

-tt E ^~~rr^~(An, NA p )f p + -^(A n , NF ± ) 

n pej p n 



(1 + A P ) 



(A n ,NAp)f p + j^(A n ,NF ± ), nej. 



(A.2) 



Inspection of equation (5) for the population inversion shows that it is at least of second 
order in the amplitudes of the active modes. This holds also for F± since the source of 
the electric field is ultimately the term in iV in the polarisation equation (4). Therefore 
the last two terms in (A.2) are respectively of third and fourth order in the amplitudes 
of the active modes. As we are interested only in a low order expansion of (A.2), we can 
drop the last term and expand the second term up to third order. 

The lowest order terms in the expansion of the nonlinear terms on the right hand side 
of (5) are 
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This suggests expanding the population inversion as 

V X -\ .;/■•//.//■• 



(A.3) 



(A.4) 



where the expansion coefficients Ajfc are such that Aj& = N k j in order to ensure that A" is 
real. Expressions for the Nj k can be obtained by substituting (A.4) into (5). Before doing 
so it is convenient to expand separately the left hand side of (5), i.e. to compute the time 
derivative of A" in terms of the time derivative of the amplitudes of the active modes: 



= y at ( f *j k + f k -f.) . 

dt 3 V df 1 ' 3 dV 3 ) 



(A.5) 



Substituting the lowest order (i.e. the linear) terms of the active modes equation (A. 2) 
into this expression we obtain 



dN 



= y: (a, + h)N jk fj k . 



(A.6) 
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Substituting (A.6) and (A.3) into (5) we obtain 

E (\- + h)N jk fjf k = - 7 E + [1 + {X ^r Xk)/2] xA 3 A k ) fj k , 



(A.7) 



from which we finally derive 
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Substituting this expression into the third term of (A. 2) and keeping in mind that the last 
term of (A. 2) can be neglected because it is of higher order we obtain the equation (21) 
for the amplitudes of the active modes. 



B Space dependent pump parameter 



If the pump parameter is space-dependent the calculations to obtain the normal forms 
become much more involved. The principle behind the derivation remains unchanged, we 
still need to separate the slow from the fast modes, but its actual realisation is a non 
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trivial affair. If the pump is flat, i.e. x does not depend on the transverse coordinates, the 
cavity modes diagonalise the linear equations: therefore the slow modes are a subset of 
the cavity modes and we can relatively easily split the laser equations into two blocks. If 
the pump parameter is space dependent this is no longer true. The eigenmodes G n (x,y) 
of the linearised laser equations (9) are no longer cavity modes, but can be expressed as 
linear superpositions of them: 

oo 

G n (x,y) = J29kA k (x,y). (B.l) 
k=i 

In order to separate the slow and fast dynamics we then need to project the equations (3-5) 
onto the modes G n using the decomposition (B.l). This is of course possible in principle, 
but rather hard to do in practice, especially since the modes G n (x,y) are usually only 
known numerically. Moreover, the laser equations are not diagonal on the basis of the 
modes G n and the resulting equations for the amplitudes of these modes contain sums 
over all the coefficients f3 n of equation (3). It is likely that a symbolic algebra package 
could be effective in splitting fast and slow variables, after which the procedure detailed 
in Appendix A can in principle be applied to obtain the normal forms. 
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